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Abstract. This paper is a survey on universal algorithms for solving the matrix Bell- 
man equations over semirings and especially tropical and idempotent semirings. However, 
original algorithms are also presented. Some applications and software implementations 
are discussed. 



1. Introduction 

Computational algorithms are constructed on the basis of certain primitive operations. 
These operations manipulate data that describe "numbers." These "numbers" are el- 
ements of a "numerical domain," i.e., a mathematical object such as the field of real 
numbers, the ring of integers, different semirings etc. 

In practice, elements of the numerical domains are replaced by their computer rep- 
resentations, i.e., by elements of certain finite models of these domains. Examples of 
models that can be conveniently used for computer representation of real numbers are 
provided by various modifications of floating point arithmetics, approximate arithmetics 
of rational numbers [30], interval arithmetics etc. The difference between mathematical 
objects ("ideal" numbers) and their finite models (computer representations) results in 
computational (e.g., rounding) errors. 

An algorithm is called universal if it is independent of a particular numerical domain 
and/or its computer representation [251 1261 [Ml 129]. A typical example of a universal 
algorithm is the computation of the scalar product (x, y) of two vectors x — (xi, . . . , x n ) 
and y = (yi, . . . , y n ) by the formula (x, y) = x±y± + • • • + x n y n . This algorithm (formula) is 
independent of a particular domain and its computer implementation, since the formula 
is well-defined for any semiring. It is clear that one algorithm can be more universal than 
another. For example, the simplest Newton-Cotes formula, the rectangular rule, provides 
the most universal algorithm for numerical integration. In particular, this formula is valid 
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also for idempotent integration (that is, over any idempotent semiring, see e.g. [201 EI]- 
Other quadrature formulas (e.g., combined trapezoid rule or the Simpson formula) are 
independent of computer arithmetics and can be used (e.g., in the iterative form) for 
computations with arbitrary accuracy. In contrast, algorithms based on Gauss- Jacob i 
formulas are designed for fixed accuracy computations: they include constants (coefficients 
and nodes of these formulas) defined with fixed accuracy. (Certainly, algorithms of this 
type can be made more universal by including procedures for computing the constants; 
however, this results in an unjustified complication of the algorithms.) 

Modern achievements in software development and mathematics make us consider nu- 
merical algorithms and their classification from a new point of view. Conventional numer- 
ical algorithms are oriented to software (or hardware) implementation based on floating 
point arithmetic and fixed accuracy. However, it is often desirable to perform computa- 
tions with variable (and arbitrary) accuracy. For this purpose, algorithms are required 
that are independent of the accuracy of computation and of the specific computer repre- 
sentation of numbers. In fact, many algorithms are independent not only of the computer 
representation of numbers, but also of concrete mathematical (algebraic) operations on 
data. In this case, operations themselves may be considered as variables. Such algorithms 
are implemented in the form of generic programs based on abstract data types that are 
defined by the user in addition to the predefined types provided by the language. The cor- 
responding program tools appeared as early as in Simula-67, but modern object-oriented 
languages (like C + +, see, e.g., [551 S3]) are more convenient for generic programming. 
Computer algebra algorithms used in such systems as Mathematica, Maple, REDUCE, 
and others are also highly universal. 

A different form of universality is featured by iterative algorithms (beginning with 
the successive approximation method) for solving differential equations (e.g., methods 
of Euler, Euler-Cauchy, Runge-Kutta, Adams, a number of important versions of the 
difference approximation method, and the like), methods for calculating elementary and 
some special functions based on the expansion in Taylor's series and continuous fractions 
(Pade approximations). These algorithms are independent of the computer representation 
of numbers. 

The concept of a generic program was introduced by many authors; for example, in [23J 
such programs were called 'program schemes.' In this paper, we discuss universal algo- 
rithms implemented in the form of generic programs and their specific features. This paper 
is closely related to [2U [251 ESI [291 EH [28], 26] , in which the concept of a universal algorithm 
was defined and software and hardware implementation of such algorithms was discussed 
in connection with problems of idempotent mathematics, see, e.g., [20 ] 1331 [38] [5T ] [52]. 
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The so-called idempotent correspondence principle, see j25j[26], linking this mathematics 
with the usual mathematics over fields, will be discussed below. In a nutshell, there 
exists a correspondence between interesting, useful, and important constructions and 
results concerning the field of real (or complex) numbers and similar constructions dealing 
with various idempotent semirings. This correspondence can be formulated in the spirit 
of the well-known N. Bohr's correspondence principle in quantum mechanics; in fact, 
the two principles are closely connected (see [2U [251 125] )• In a sense, the traditional 
mathematics over numerical fields can be treated as a 'quantum' theory, whereas the 
idempotent mathematics can be treated as a 'classical' shadow (or counterpart) of the 
traditional one. It is important that the idempotent correspondence principle is valid for 
algorithms, computer programs and hardware units. 

In quantum mechanics the superposition principle means that the Schrodinger equa- 
tion (which is basic for the theory) is linear. Similarly in idempotent mathematics the 
(idempotent) superposition principle (formulated by V. P. Maslov) means that some im- 
portant and basic problems and equations that are nonlinear in the usual sense (e.g., the 
Hamilton- Jacobi equation, which is basic for classical mechanics and appears in many 
optimization problems, or the Bellman equation and its versions and generalizations) can 
be treated as linear over appropriate idempotent semirings, see [36l [37] . 

Note that numerical algorithms for infinite dimensional linear problems over idem- 
potent semirings (e.g., idempotent integration, integral operators and transformations, 
the Hamilton- Jacobi and generalized Bellman equations) deal with the corresponding 
finite-dimensional approximations. Thus idempotent linear algebra is the basis of the 
idempotent numerical analysis and, in particular, the discrete optimization theory. 

B. A. Carre [TJE] ( see also [112 EH E]) used the idempotent linear algebra to show that 
different optimization problems for finite graphs can be formulated in a unified manner 
and reduced to solving Bellman equations, i.e., systems of linear algebraic equations over 
idempotent semirings. He also generalized principal algorithms of computational linear 
algebra to the idempotent case and showed that some of these coincide with algorithms 
independently developed for solution of optimization problems. For example, Bellman's 
method of solving the shortest path problem corresponds to a version of Jacobi's method 
for solving a system of linear equations, whereas Ford's algorithm corresponds to a version 
of Gauss-Seidel's method. We briefly discuss Bellman equations and the corresponding 
optimization problems on graphs, and use the ideas of Carre to obtain new universal 
algorithms. We stress that these well-known results can be interpreted as a manifestation 
of the idempotent superposition principle. 
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Note that many algorithms for solving the matrix Bellman equation could be found 
in^H3ISl|Tni|T51inil2Sll2SllSlllSll5D]. More general problems of linear algebra over the 
max-plus algebra are examined, e.g. in [5]. 

We also briefly discuss interval analysis over idempotent and positive semirings. Idem- 
potent interval analysis appears in (32J [331 HE], where it is applied to the Bellman matrix 
equation. Many different problems coming from the idempotent linear algebra, have been 
considered since then, see e.g. (HI [321 E333 EU HI]- It is important to observe that intervals 
over an idempotent semiring form a new idempotent semiring. Hence universal algorithms 
can be applied to elements of this new semiring and generate interval extensions of the 
initial algorithms. 

This paper is about software implementations of universal algorithms for solving the 
matrix Bellman equations over semirings. In Section [2] we present an introduction to math- 
ematics of semirings and especially to the tropical (idempotent) mathematics, i.e., the area 
of mathematics working with idempotent semirings (that is, semirings with idempotent 
addition). In Section [3] we present a number of well-known and new universal algorithms 
of linear algebra over semirings, related to discrete matrix Bellman equation and algebraic 
path problem. These algorithms are closely related to their linear-algebraic prototypes 
described, e.g., in the celebrated book of Golub and Van Loan [14] which serves as the 
main source of such prototypes. Following the style of [H] we present them in MATLAB 
code. The perspectives and experience of their implementation are also discussed. 

2. Mathematics of semirings 

2.1. Basic definitions. A broad class of universal algorithms is related to the concept of 
a semiring. We recall here the definition of a semiring (see, e.g., [13] )• Consider a semiring, 
i.e., a set S endowed with two associative operations: addition © and multiplication 
such that addition is commutative, multiplication distributes over addition from either 
side, (resp., 1) is the neutral element of addition (resp., multiplication), OOx = xQO = 
for all x E S, and 0^1. Let the semiring S be partially ordered by a relation ^ such that 
is the least element and the inequality x <y implies that xQ)z^yQ)z,xQz^y(Dz, 
and z x ^ z y for all x, y, z G S; in this case the semiring S is called positive (see, 
e.g., [13]). 

A semiring S is called a semifield if every nonzero element is invertible. 

A semiring S is called idempotent if X X — X for all x £ S. In this case the addition © 
defines a canonical partial order ^ on the semiring S by the rule: x ^ y iff x © y = y. It 
is easy to prove that any idempotent semiring is positive with respect to this order. Note 
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also that x © y = sup{x,y} with respect to the canonical order. In the sequel, we shall 
assume that all idempotent semirings are ordered by the canonical partial order relation. 

We shall say that a positive (e.g., idempotent) semiring S is complete if it is complete 
as an ordered set. This means that for every subset T C S there exist elements supT G S 
and inf T G S. 

The most well-known and important examples of positive semirings are "numerical" 
semirings consisting of (a subset of) real numbers and ordered by the usual linear order 
< on R: the semiring R + with the usual operations © = +,© = ■ and neutral elements 
= 0, 1 = 1, the semiring R max = R U {— oo} with the operations © = max, = + 
and neutral elements = — oo, 1 = 0, the semiring R max = R ma x U {°°}) where x ^ oo, 
x © oo = oo for all x, x oo = oo © x = oo if x ^ 0, and oo = oo 0, and the 
semiring S^l min = [a, b], where — oo < a < b < +oo, with the operations © = max, 
= min and neutral elements = a, 1 = b. The semirings R max , Rmax, and m - m = 
[a,b] are idempotent. The semirings R max , 'S'maxmin; ^+ = ^+ U{°°} are complete as 
ordered sets. Remind that every partially ordered set can be imbedded to its completion 
(a minimal complete set containing the initial one). The semiring R rmn = R|J{oo} with 
operations © = min and = + and neutral elements = oo, 1 = is isomorphic to 



The semiring R max is also called the max-plus algebra. The semifields R max and R m i n 
are called tropical algebras. The term "tropical" initially appeared in [UJ for a discrete 
version of the max-plus algebra as a suggestion of Ch. Choffrut, see also [TSJ EHJ 152] - 

Many mathematical constructions, notions, and results over the fields of real and com- 
plex numbers have nontrivial analogs over idempotent semirings. Idempotent semirings 
have become recently the object of investigation of new branches of mathematics, idem- 
potent mathematics and tropical geometry, see, e.g. [21 [TUl [2U EH EU 152] . 

Denote by Mai mn (S) a set of all matrices A = (ay) with m rows and n columns whose 
coefficients belong to a semiring S. The sum A © B of matrices A,Be M&t mn (S) and 
the product AB of matrices A e Mat/ m (S') and B e Mat mn (5') are defined according to 
the usual rules of linear algebra: A © B = (oy © fey) G Mat mn (S') and 



where A 6 Mat im (S') and B G Mat mn (5'). Note that we write AB instead of A B. 

If the semiring S is positive, then the set 
Mat mn (S) is ordered by the relation A = (a^) H B = (pij) iff -< bij in S for all 
1 < i < m, 1 < j < n. 
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The matrix multiplication is consistent with the order ^ in the following sense: if 
A, A' G Mat lm (S), B, B' G Mat mn (5) and A ^ A', B H B', then AB ^ A'B' in M&t ln (S). 
The set Mat nn (S') of square {n x n) matrices over a [positive, idempotent] semiring S 
forms a [positive, idempotent] semi-ring with a zero element O = (o^), where oy = 0, 
1 < 2, j ; < n, and a unit element / = (Sij), where 8{j = 1 if i = j and 5ij = otherwise. 

The set Mat nn is an example of a noncommutative semiring if n > 1. 

2.2. Closure operations. Let a positive semiring S be endowed with a partial unary 
closure operation * such that a implies a* -< b* and a* = 1 © (a* a) = 1 © (a a*) 
on its domain of definition. In particular, 0* = 1 by definition. 

These axioms imply that a* = 1 © a © a 2 © • • • © (a* a") if n > 1. Thus x* can be 
considered as a 'regularized sum' of the series a* = l©a©a 2 ©. . . . In a positive semiring, 
provided that it is closed under taking bounded ordered sups and the operations © and 
distribute over such sups we can define 



(1) a* := sup 1 © a © ... © a k , 

k>0 

if the sequence on the r.h.s. is bounded. In this case a* is the least solution of x — ax©l 
and x = xaffil, and more generally a*b is the the least solution of x = ax®b and x = xa®b. 
So if S is complete, then the closure operation is well-defined for every element x G S. 
In the case of idempotent addition ([T|) becomes particularly nice: 

(2) a* = £fV = sup a*. 

So ^° 

In numerical semirings the operation * is usually very easy to implement: x* = (1 — x)^ 1 
if x < 1 in R + , or R + and x* = oo if x > 1 in R + ; x* = 1 if x ^ 1 in R max and R max , 
x* = oo if x y 1 in R max , x* = 1 for all x in S^ia^mm- I* 1 au other cases x* is undefined. 

The closure operation in matrix semirings over a positive semiring S can be defined 
inductively: A* = (an)* = (a* x ) in Matn(S) and for any integer n > 1 and any matrix 

'A n A 12 
A21 A22 

where An G Mat fcfe (5'), A i2 G Mat fcn _ Jt (S'), A 2 i G Mat n - kk (S), A 22 G Mat^-fen-A^S 1 ), 
1 < k < n, by defintion, 

/ A* n © A n A 12 D*A 21 A* n A* n A 12 D* 1 

(3) A* = 

\ D*A 21 A* U D* 

where D = A22 © ^42i^4*i^i2- bt can be proved that this definition of A* implies that the 
equalities A* = A* A® I = AA* ©/ are satisfied, and thus A* is a 'regularized sum' of the 
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series I © A © A 2 © ... . Moreover, in the case when A* is defined as the least solution of 
A* = A* A © I and A* = AA* © I, it can be shown that it satisfies ©. 

Note that this recurrence relation coincides with the formulas of escalator method of 
matrix inversion in the traditional linear algebra over the field of real or complex numbers, 
up to the algebraic operations used. Hence this algorithm of matrix closure requires a 
polynomial number of operations in n, see below for more details. 

Let S be a positive semiring. The matrix (or discrete stationary) Bellman equation has 
the form 

(4) X = AX © B, 

where A G Mat nn (5'), I,B £ Mat ns (S'), and the matrix X is unknown. Let A* be the 
closure of the matrix A. It follows from the identity A* = A* A © / that the matrix A*B 
satisfies this equation. As in the scalar case, it can be shown that for positive semirings 
under reasonable distributivity assumptions, if A* is defined as in (pQ) then A*B is the 
least in the set of solutions to equation (TJJ with respect to the partial order in Mat ns (5'). 
Recall that in the idempotent case 

(5) A* = 0A ! = sup A 1 . 

So ^° 

Consider also the case when A = (ay) is n x n strictly upper-triangular (such that 
Oij = for i > j) , or n x n strictly lower-triangular (such that = for i < j). In this 
case A n = O, the all- zeros matrix, and it can be shown by iterating X = AX © I that 
this equation has a unique solution, namely 

(6) A* = I@A@...®A n -\ 

Curiously enough, formula (jH]) works more generally in the case of numerical idempotent 
semirings: in fact, the series ([5]) converges there if and only if it can be truncated to ([6]). 
This is closely related to the principal path interpretation of A* explained in the next 
subsection. 

2.3. Weighted directed graphs and matrices over semirings. Suppose that S is a 
semiring with zero and unity 1. It is well-known that any square matrix A = (dy) G 
Mat nn (S') specifies a weighted directed graph. This geometrical construction includes 
three kinds of objects: the set X of n elements xi,...,x n called nodes, the set T of 
all ordered pairs (xi,Xj) such that ^ called arcs, and the mapping A: Y — > S 
such that A(xi The elements of the semiring S are called weights of the 

arcs. Conversely, any given weighted directed graph with n nodes specifies a unique matrix 
A G Mat nn (S). 
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Figure 1. 




This definition allows for some pairs of nodes to be disconnected if the corresponding 
element of the matrix A is and for some channels to be "loops" with coincident ends if 
the matrix A has nonzero diagonal elements. 

Recall that a sequence of nodes of the form 

V = {yo,Vi, ■■■,Vk) 

with k > and j/i+i) G T, i — 0, . . . , k — 1, is called a path of length connecting 
y with 2/&. Denote the set of all such paths by Pkiuo^Uk)- The weight A{p) of a path 
p G Pk(yo, Vk) is defined to be the product of weights of arcs connecting consecutive nodes 
of the path: 

A(p) = A(y , yi) © • ■ ■ A(y k _ u y k ). 

By definition, for a 'path' p G Po(xi,Xj) of length k = the weight is 1 if i — j and 
otherwise. 

For each matrix A G Mat n „(S') define A = E = (5^-) (where 5^ = 1 if i — j and 
= otherwise) and A k = AA k ~ 1 , k > 1. Let a'y' be the (z, j)th element of the matrix 
A k . It is easily checked that 

a ij = a «o*i © ' ' ' © a *fc-i*fc- 

«0=«, ik=j 
l<h,...,i k _ 1 <n 

Thus is the supremum of the set of weights corresponding to all paths of length k 
connecting the node Xi = Xi with xi k = Xj. 

Let A* be defined as in ((SJ). Denote the elements of the matrix A* by a|-, i, j = 1, . . . , n; 
then 

0<fc<oopgP fc (x i ,a: i ) 

The closure matrix A* solves the well-known algebraic path problem, which is formulated 
as follows: for each pair (x{, Xj) calculate the supremum of weights of all paths (of arbitrary 
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length) connecting node X{ with node x r The closure operation in matrix semirings has 
been studied extensively (see, e.g., |2j [7J [101 [131 EH CEO [2Ql [33] and references therein). 

Example 1 (The shortest path problem). Let S = R m in, so the weights are real numbers. 
In this case 

A(p) = A(y Q ,yx) + A(y 1 ,y 2 ) H h A(y k ^,y k ). 

If the element specifies the length of the arc (xi,Xj) in some metric, then a*j is the 
length of the shortest path connecting x-i with x r 

Example 2 (The maximal path width problem). Let S = R U {0, 1} with © = max, 
= min. Then 

a*- = max A(p) 1 

PG U Pk{Xi,Xj) 
k>l 

A(p) = mm(A(y , y{), . . . , A(y k -i, J/*))- 

If the element a^- specifies the "width" of the arc (xi,Xj), then the width of a path p 
is defined as the minimal width of its constituting arcs and the element a*- gives the 
supremum of possible widths of all paths connecting Xi with Xj. 

Example 3 (A simple dynamic programming problem). Let S = R max and suppose 
a>ij gives the profit corresponding to the transition from Xi to Xj. Define the vector 
B = (b{) G Mat nl (R max ) whose element 6j gives the terminal profit corresponding to 
exiting from the graph through the node X{. Of course, negative profits (or, rather, 
losses) are allowed. Let m be the total profit corresponding to a path p e Pk{%i, i-e. 

m = A{p) + by 

Then it is easy to check that the supremum of profits that can be achieved on paths 
of length k beginning at the node is equal to (A k B)i and the supremum of profits 
achievable without a restriction on the length of a path equals (A*B){. 

Example 4 (The matrix inversion problem). Note that in the formulas of this section 
we are using distributivity of the multiplication © with respect to the addition © but 
do not use the idempotency axiom. Thus the algebraic path problem can be posed for a 
nonidempotent semiring S as well (see, e.g., fS]). For instance, if S = R, then 

A* = I + A + A 2 + ■ ■ ■ = (/ - A)' 1 . 

If \\A\\ > 1 but the matrix I — A is invertible, then this expression defines a regularized 
sum of the divergent matrix power series Yli>o^ - 
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We emphasize that this connection between the matrix closure operation and solutions 
to the Bellman equation gives rise to a number of different algorithms for numerical 
calculation of the matrix closure. All these algorithms are adaptations of the well-known 
algorithms of the traditional computational linear algebra, such as the Gauss- Jordan 
elimination, various iterative and escalator schemes, etc. This is a special case of the 
idempotent superposition principle (see below). 

In fact, the theory of the discrete stationary Bellman equation can be developed using 
the identity A* = AA* © / as an additional axiom without any substantial interpretation 
(the so-called closed semirings, see, e.g., [131 1221 E]). Such theory can be based on the 
following identities, true both for the case of idempotent semirings with path interpre- 
tation, and the real numbers with conventional arithmetic (assumed that A and B have 
appropriate sizes): 

(A® BY = (A*B)*A*, 

(7) 

(AB)*A = A(BA)* . 

2.4. Interval analysis over positive semirings. Traditional interval analysis is a non- 
trivial and popular mathematical area, see, e.g., [U [121 I2U Ell H2]. An "idempotent" 
version of interval analysis (and moreover interval analysis over positive semirings) ap- 
peared in [321 [331 EE]- Rather many publications on the subject appeared later, see, 
e.g., [121 ESI HQl HE]- Interval analysis over the positive semiring R + was discussed 
in [4]. 

Let a set S be partially ordered by a relation -<. A closed interval in S is a subset of 
the form x= [x,x] = {a; 6 5 | a; ^x}, where the elements x ■< x are called lower 
and upper bounds of the interval x. The order -< induces a partial ordering on the set of 
all closed intervals in S: x -< y iff x -< y and x -< y. 

A weak interval extension I(S) of a positive semiring S is the set of all closed intervals 
in S endowed with operations © and defined asx©y = [xffiy, x©y], x0y = 
[x y, x y] and a partial order induced by the order in S. The closure operation in 
I(S) is defined by x* = [x*,x*]. There are some other interval extensions (including the 
so-called strong interval extension [33]) but the weak extension is more convenient. 

The extension I(S) is positive; I(S) is idempotent if S is an idempotent semiring. A 
universal algorithm over S can be applied to I(S) and we shall get an interval version of 
the initial algorithm. Usually both versions have the same complexity. For the discrete 
stationary Bellman equation and the corresponding optimization problems on graphs, 
interval analysis was examined in [321 [33] in details. Other problems of idempotent linear 
algebra were examined in P H21 ESI SD1 BE] ■ 
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Idempotent mathematics appears to be remarkably simpler than its traditional analog. 
For example, in traditional interval arithmetic, multiplication of intervals is not distribu- 
tive with respect to addition of intervals, whereas in idempotent interval arithmetic this 
distributivity is preserved. Moreover, in traditional interval analysis the set of all square 
interval matrices of a given order does not form even a semigroup with respect to matrix 
multiplication: this operation is not associative since distributivity is lost in the traditional 
interval arithmetic. On the contrary, in the idempotent (and positive) case associativity 
is preserved. Finally, in traditional interval analysis some problems of linear algebra, such 
as solution of a linear system of interval equations, can be very difficult (more precisely, 
they are ./VP-hard, see [2T] and references therein). It was noticed in [521 155] that in the 
idempotent case solving an interval linear system requires a polynomial number of opera- 
tions (similarly to the usual Gauss elimination algorithm). Two properties that make the 
idempotent interval arithmetic so simple are monotonicity of arithmetic operations and 
positivity of all elements of an idempotent semiring. 

Interval estimates in idempotent mathematics are usually exact. In the traditional 
theory such estimates tend to be overly pessimistic. 

2.5. Idempotent correspondence principle. There is a nontrivial analogy between 
mathematics of semirings and quantum mechanics. For example, the field of real num- 
bers can be treated as a "quantum object" with respect to idempotent semirings. So 
idempotent semirings can be treated as "classical" or "semi-classical" objects with re- 
spect to the field of real numbers. 

Let R be the field of real numbers and R + the subset of all non-negative numbers. 
Consider the following change of variables: 

u i — y w — h In u, 

where u G R+ \ {0}, h > 0; thus u = e w l } \ «j£R. Denote by the additional element 
— oo and by S the extended real line RU{0}. The above change of variables has a natural 
extension D h to the whole S by D h (0) = 0; also, we denote D h {l) = = 1. 

Denote by Sh the set S equipped with the two operations Qh (generalized addition) 
and Qh (generalized multiplication) such that Dh is a homomorphism of {R + ,+,-} to 
{S, Q h , Q h }. This means that D h {u x +u 2 ) = D h (ui)® h D h (u 2 ) and D h (u v u 2 ) = D h {ui)Q h 
Dh(u 2 ), i.e., w\ Qh w 2 = wi + w 2 and W\ Qh w 2 = hln(e Wl ^ h + e W2 ^ h ). It is easy to prove 
that W\ Qh w 2 — » max{w 1; w 2 } as h — > 0. 

Algebraic structures in R + and Sh are isomorphic; therefore we have obtained R max as 
a result of a deformation of the structure in R + . We stress the obvious analogy with the 
quantization procedure, where h is the analog of the Planck constant. In these terms, R + 
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(or R) plays the part of a "quantum object" while R max acts as a "classical" or "semi- 
classical" object that arises as the result of a dequantization of this quantum object. In 
the case of R m i n , the corresponding dequantization procedure is generated by the change 
of variables u ^ w = — hhau. 

There is a natural transition from the field of real numbers or complex numbers to the 
idempotent semiring R max (or R m i n ). This is a composition of the mapping x h-> \x\ and 
the deformation described above. 

In general an idempotent dequantization is a transition from a basic field to an idempo- 
tent semiring in mathematical concepts, constructions and results, see [2H [26] for details. 
Idempotent dequantization suggests the following formulation of the idempotent corre- 
spondence principle: 

There exists a heuristic correspondence between interesting, useful, and im- 
portant constructions and results over the field of real (or complex) num- 
bers and similar constructions and results over idempotent semirings in the 
spirit of N. Bohr's correspondence principle in quantum mechanics. 



TRADITIONAL 
MATHEMATICS 



Idempotent Correspondence Principle 



I 



IDEMPOTENT 
MATHEMATICS 




QUANTUM 
MECHANICS 


N. Bohr's Correspondence Principle 


CLASSICAL 
MECHANICS 





Thus idempotent mathematics can be treated as a "classical shadow (or counterpart)" 
of the traditional Mathematics over fields. A systematic application of this correspondence 
principle leads to a variety of theoretical and applied results, see, e.g. 
EH [52]. Relations to quantum physics are discussed in detail, e.g., in 

In this paper we aim to develop a practical systematic application of the correspondence 
principle to the algorithms of linear algebra and discrete mathematics. For the remainder 
of this subsection let us focus on an idea how the idempotent correspondence principle 
may lead to a unifying approach to hardware design. (See [SUES] for more information.) 

The most important and standard numerical algorithms have many hardware realiza- 
tions in the form of technical devices or special processors. These devices often can be used 
as prototypes for new hardware units resulting from mere substitution of the usual arith- 
metic operations by their semiring analogs (and additional tools for generating neutral 
elements and 1). Of course, the case of numerical semirings consisting of real numbers 
(maybe except neutral elements) and semirings of numerical intervals is the most simple 
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and natural. Note that for semifields (including R max and R m i n ) the operation of division 
is also defined. 

Good and efficient technical ideas and decisions can be taken from prototypes to new 
hardware units. Thus the correspondence principle generates a regular heuristic method 
for hardware design. Note that to get a patent it is necessary to present the so-called 
'invention formula', that is to indicate a prototype for the suggested device and the 
difference between these devices. 

Consider (as a typical example) the most popular and important algorithm of computing 
the scalar product of two vectors: 

(8) (x, y) = xiyi + x 2 y2 H h x n y n . 

The universal version of jS} for any semiring A is obvious: 

(9) (x, y) = (xi yi) © (x 2 y 2 ) © • • • © (x n © y n )- 
In the case A = R max this formula turns into the following one: 

(10) (x, y) = max{xi + y 1} x 2 + y 2 ,--- , x n + y n }- 

This calculation is standard for many optimization algorithms, so it is useful to con- 
struct a hardware unit for computing (1101) . There are many different devices (and patents) 
for computing (jSJ) and every such device can be used as a prototype to construct a new 
device for computing (fit)]) and even ([9]). Many processors for matrix multiplication and 
for other algorithms of linear algebra are based on computing scalar products and on the 
corresponding "elementary" devices. Using modern technologies it is possible to construct 
cheap special-purpose multi-processor chips and systolic arrays of elementary processors 
implementing universal algorithms. See, e.g., [3U |2SJ IB] where the systolic arrays and 
parallel computing issues are discussed for the algebraic path problem. In particular, 
there is a systolic array of n(n + 1) elementary processors which performs computations 
of the Gauss- Jordan elimination algorithm and can solve the algebraic path problem 
within 5n — 2 time steps. 

3. Some universal algorithms of linear algebra 

In this section we discuss universal algorithms computing A* and A*B. We start with 
the basic escalator and Gauss- Jordan elimination techniques in Subsect. 13. II and continue 
with its specification to the case of Toeplitz systems in Subsect. 13. 21 The universal 
LDM decomposition of Bellman equations is explained in Subsect. I3.3[ followed by its 
adaptations to symmetric and band matrices in Subsect. 13.41 The iteration schemes are 
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discussed in Subsect. 13.51 In the final Subsect. 13.61 we discuss the implementations of 
universal algorithms. 

Algorithms themselves will be described in a language of Matlab, following the tradition 
of Golub and van Loan [14|. This is done for two purposes: 1) to simplify the comparison 
of the algorithms with their prototypes taken mostly from [13], 2) since the language of 
Matlab is designed for matrix computations. We will not formally describe the rules of 
our Matlab-derived language, preferring just to outline the following important features: 

1. Our basic arithmetic operations are a © b, a b and a*. 

2. The vectorization of these operations follows the rules of Matlab. 

3. We use basic keywords of Matlab like 'for', 'while', 'if and 'end', similar to other 
programming languages like C + + or Java. 

Let us give some examples of universal matrix computations in our language: 
Example 1. v(l : j) = a* a(l : j, k) means that the result of (scalar) multiplication of 
the first j components of the kth column of A by the closure of a is assigned to the first 
j components of v. 

Example 2. a(i,j) = a(i,j) © a(i, 1 : n) a(l : n,j) means that we add (©) to the entry 
ay of A the result of the (universal) scalar multiplication of the zth row with the jth 
column of A (assumed that A is n x n). 

Example 3. a(l : n, k) b(l, 1 : m) means the outer product of the kth column of A with 
the Zth row of B. The entries of resulting matrix C = (cy) equal = by, for all 
i — 1, . . . , n and j = 1, . . . , m. 

Example 4- x(l : n) y(n : — 1 : 1) is the scalar product of vector x with vector 
y whose components are taken in the reverse order: the proper algebraic expression is 

©i=l x i © Vn+l—i- 

Example 5. The following cycle yields the same result as in the previous example: s = 
for i = 1 : n 

s = s © x(i) x(n + 1 — i) 



3.1. Escalator scheme and Gauss- Jordan elimination. We first analyse the basic 
escalator method, based on the definition of matrix closures Let A be a square matrix. 
Closures of its main submatrices A^ can be found inductively, starting from A\ = (an)*, 
the closure of the first diagonal entry. Generally we represent Ak+i as 



end 
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assuming that we have found the closure of A}.. In this representation, g k and hk are 
columns with k entries and a k+ i is a scalar. We also represent A* k+1 as 



Al+i — l T 



wT U k+ i 



Using (jHJ) we obtain that 



(11) 



u k+i = {hlA* k g k © a k+1 ] 
v k = A* k g k u k+1 , 

.t _ „. r T A* 



w k = u k+1 h k A* k , 
U k = Alg k u k+l h T k Al © A*. 
An algorithm based on (ITT]) can be written as follows. 

Algorithm 1. Escalator method for computing A* 

Input: an n x n matrix A with entries a(i,j), 
also used to store the final result 

and the intermediate results of the computation process. 

a(l,l) = (a(l,l))* 
for i — 1 : n — 1 

Ag = a(l : i, 1 : z) a(l : z, z + 1) 
= a(i + 1, 1 : i) a(l : i, 1 : i) 
a(i + 1, i + 1) = a(i + 1, i + 1) © a(i + 1, 1 : i) © : z, 1) 
a(z' + l,z + l) = (a(z + l,z' + l))* 
a(l :i,i + l) = a(i + l,i + l)0As 
a(i + 1, 1 : i) = a(i + 1, % + 1) hA 
a(l : z, 1 : z) = a(l : z, 1 : z) © a(z + 1, z + 1) M 
end 

In full analogy with its linear algebraic prototype, the algorithm requires n 3 + 0(n 2 ) 
operations of addition ©, n 3 + 0(n 2 ) operations of multiplication ©, and n operations of 
taking algebraic closure. The linear-algebraic prototype of the method written above is 
also called the bordering method in the literature [TTj . 

Alternatively, we can obtain a solution of X = AX © B as a result of elimination 
process, whose informal explanation is given below. If A* is defined as © i>0 A 1 (including 
the scalar case), then A*B is the least solution of X = AX © B for all A and B of 
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appropriate sizes. In this case, the solution found by the elimination process given below 
coincides with A*B. 

For matrix A = (a^-) and column vectors x = (xi) and b = (hi) (restricting w.l.o.g. to 
the column vectors), the Bellman equation x = Ax © b can be written as 



(12) 



•?'2 



a>i2 



0,21 a 22 



\ a nl a n2 



a h 

&2n 



\ 



'J 



x 2 



(l o ... o\ 

1 ... 



\0 



7 



\b n j 



After expressing x\ in terms of x 2 , ■ ■ ■ ,x n from the first equation and substituting this 
expression for X\ in all other equations from the second to the nth we obtain 



x 2 



[an a 12 



/o 

a 2 2 © (a 2 i(an)*ai 2 ) 



(an) air 



\ 



a 2n © (a 2 i(an)*ai n ) 



(13) 



\x n ) \0 a n2 © (a n i(au)*ai 2 ) 

/ (an)* 
«2i (an)* 1 



/xi\ 



x 2 



a nn © (a n i(aii)*ai n )y \x n J 
0\ fh\ 







\bn) 



\a„i(an)* ... 1/ 

Note that nontrivial entries in both matrices occupy complementary places, so during 
computations both matrices can be stored in the same square array C^. Denote its 
elements by where k is the number of eliminated variables. After I — 1 eliminations 
we have 



(14) 





= (<r i, )x 










c (0 
c il 






;l- 


1,/ + 1,.. 


. ,n 




= 4r i, ®<=r i) (4 , - l) ; 


*<t 


1) 

j 








= 1,...,/-1,(+1,.. 


. ,n 










= (4 , - i, )*cr ,) , < = 


1,.. 




1,1 + 1,.. 


.,n 



After n eliminations we K et x = C^b. Taking as 6 any vector with one coordinate equal 
to 1 and the rest equal to 0, we obtain C^ n ' = A*. We write out the following algorithm 
based on recursion (JHJ). 
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Algorithm 2. Gauss- Jordan elimination for computing A*. 
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Input: annxn matrix A with entries a(i,j), 

also used to store the final result 

and intermediate results of the computation process. 

for i — 1 : n 

a(i,i) = (a(i,i))* 
for k = 1 : n 
if k 

a(k, i) = a(k, i) a(i, i) 

end 

end 

for k — 1 : n 
for j = 1 : n 

iikj^iSzj^i 

a (k, j) = a(k, j) © a(k, i) a(i, j) 

end 

end 

for j — 1 : n 
if j i 

a(hj) = a(ij) 

end 

end 

end 

Remark 1. Algorithm [2] can be regarded as a "universal Floyd-Warshall algorithm" 
generalizing the well-known algorithms of Warshall and Floyd for computing the transitive 
closure of a graph and all optimal paths on a graph. See, e.g., [15] for the description of 
these classical methods of discrete mathematics. In turn, these methods can be regarded 
as specifications of Algorithm [2] to the cases of max-plus and Boolean semiring. 

Remark 2. Algorithm[2]is also close to Yershov's "refilling" method for inverting matrices 
and solving systems Ax = b in the classical linear algebra, see [11] Chapter 2 for details. 

3.2. Toeplitz systems. We start by considering the escalator method for finding the 
solution x = A*b to x = Ax © b, where x and b are column vectors. Firstly, we have 
x^ = A\b\. Let x^ be the vector found after (k — 1) steps, and let us write 
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X 



(fc+1) 




(15) 



Using pip we obtain that 

= u k+1 (hlx {k) © b k+ i), 
z = x {k) © A* k g k x k+1 . 

We have to compute A* k g k . In general, we would have to use Algorithm [U Next we 
show that this calculation can be done very efficiently when A is symmetric Toeplitz. 

Formally, a matrix A G Mat„„(<S) is called Toeplitzif there exist scalars r_ n +i, . . . , ro, . . . , r n _i 



such that Aj 



1 3-1 



for all i and j. Informally, Toeplitz matrices are such that their entries 



are constant along any line parallel to the main diagonal (and along the main diagonal 
itself). For example, 



/ 



A 



r n 
r_i r 
r_ 2 r_i 



r 2 
r 



'3 

''i 



yr_ 3 r_ 2 r_i r y 

is Toeplitz. Such matrices are not necessarily symmetric. However, they are always 
persymmetric, that is, symmetric with respect to the inverse diagonal. This property is 
algebraically expressed as A = E n A T E n , where E n = [e n , . . . ,e%]. By we denote the 
column whose iih entry is 1 and other entries are O. The property E% = I n (where I n 
is the n x n identity matrix) implies that the product of two persymmetric matrices is 
persymmetric. Hence any degree of a persymmetric matrix is persymmetric, and so is the 
closure of a persymmetric matrix. Thus, if A is persymmetric, then 



(16) 



E n A* = (A*) T E n 



Further we deal only with symmetric Toeplitz matrices. Consider the equation y 



and T n is defined by the scalars r , r 1; . . . , r n _i so 



T n y © r^ n \ where = (r 1; . . . , r r 

that Tij = ru-i\ for all i and j. This is a generalization of the Yule- Walker problem [Tj 
Assume that we have obtained the least solution y^ to the system y = T k y © for 
some k such that 1 < k < n — 1, where T k is the main k x k submatrix of T n . We write 
T k+ i as 

T k E k T<to\ 
r^ T E k r 

We also write y( fc+1 ) and r^ k+1 ^ as 



Ti 



k+l 
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y(k + l) = ( Z ] r (fc+l) = | ^ 

Using ffl3|) . ffToD and the identity T fc M fc ) = y^ k \ we obtain that 

a k = (r © r^yWj'CrW^yW © r fc+1 ), 
z = E k y^a k ®y { - k \ 

Denote = ro © r^ T y^ k \ The following argument shows that /3& can be found recur- 
sively if (/SjjLi) -1 exists. 



& = r © [r (fc - 1)T r fe ] 



r © r^V^ © (r^^jW*^© 



(17) r fe K_! = © (f3* k _ 1 )- 1 a|_ x . 

Existence of (/3k_i) _1 is not universal, and this will make us write two versions of our 
algorithm, the first one involving (I17p . and the second one not involving it. We will write 
these two versions in one program and mark the expressions which refer only to the first 
version or to the second one by the MATLAB-style comments %1 and %2, respectively. 
Collecting the expressions for f3 k , a k and z we obtain the following recursive expression 
for t/W; 



Pk = r ®r^ T y( k \ %2 
/3 fc = /3 fc _ 1 ©(^_ 1 )- 1 0aL 1 , %1 
(18) a k = {hT®{r^ T E k y^®r k+1 ), 

yi k + i )= ^EkV {k) a k ®y^y 

Recursive expression (TT81) is a generalized version of the Durbin method for the Yule- 
Walker problem, see [H] Algorithm 4.7.1 for a prototype. 

Algorithm 3. The Yule- Walker problem for the Bellman equations with symmetric Toeplitz 
matrix. 
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Input: tq: scalar, 
r: n — 1 X 1 vector; 



y(l) = r*0r(l) 
P = r %1 
a = rg © r(l) 



for k = 1 : n — 1 



P = r © r(l : A;) © y(l : fc) %2 

/3 = /3 © © a 2 %1 

a = /?* (r(A; : -1 : 1) : k) © r(A; + 1)) 
z(l : fc) = y(l : k) © : -1 : 1) a 
y(l : k) = z(l : k) 
y(k + 1) = a 



end 

Output: vector y. 

In the general case, the algorithm requires 3/2n 2 + 0(n) operations © and each, and 
just n 2 + 0{n) of © and © if inversions of algebraic closures are allowed (as usual, just n 
such closures are required in both cases). 

Now we consider the problem of finding x^ = T*b^ where T n is as above and b^ n ' = 
(pi, . . . , b n ) is arbitrary. We also introduce the column vectors y^ k > which solve the Yule- 
Walker problem: y^ = T^r^ k \ The main idea is to find the expression for a;( fc+1 ) = 
T k+1 b^ k+v> involving x^> and y^ k \ We write and M fc+1 ) as 



Making use of the persymmetry of T£ and of the identities T£b k = and T£r k = y( k \ 
we specialize expressions (fT5j) and obtain that 



The coefficient r © r^ k ' T y^ = f3 k can be expressed again as = k -i © (P^-i) 1 © 
(afc-i) 2 , if the closure (p k -i)* is invertible. Using this we obtain the following recursive 




ilk = (ro © r^ T y^Y (r^ T E k x^ © b k+l ) 
v = E k y {k) fj, k © x {k) . 



expression: 
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(19) 



x 



p k = ro(Br (k)T y (k)^ %2 

P k = /3 fc _i © (^.i) -1 
fi k = (51 (r^ T E k x^ © b k+1 ) 
E k y^Hk © x^ 



(fc+i) 



Expressions ([18]) and ([19]) yield the following generalized version of the Levinson al- 
gorithm for solving linear symmetric Toeplitz systems, see [2] Algorithm 4.7.2 for a 
prototype: 



Algorithm 4. Bellman system with symmetric Toeplitz matrix 



Input: ro: scalar, 

r: lxii-1 row vector; 

b: n x 1 column vector. 

y (l) =r o *0r(l); x(l) = r* 6(1); 
/3 = r %1 
a = 7-q © r(l) 
for = 1 : n — 1 

/3 = r ©r(l : fc) © y(l : fc) %2 
/3 = /3 © (/3*)- 1 © a 2 %1 



H = p*® (r(k 
v(l : fc) = x(l 
x(l : fc) = v(l 

x(k + 1) = fJL 

if fc < n — 1 

a = /3* © ( r (fc 

z(l:A;)=y(l 
: fc) = z(l 
y(k + 1) = a 
end 
end 



-1 : 1)0 x(l : k)®b(k + l)) 

k) © y(k : -1 : 1) 0/i 

k) 



-1 : 1)0 2/(1 : fc) © r(fc + 1)) 

k)®y(k : -1 : 1) ©a 

k) 



Output: vector x. 
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In the general case, the algorithm requires 5/2n 2 + 0(n) operations © and © each, and 
just 2n 2 + 0{n) of © and © if inversions of algebraic closures are allowed (as usual, just 
n such closures are required in both cases). 

3.3. LDM decomposition. Factorization of a matrix into the product A = LDM, 
where L and M are lower and upper triangular matrices with a unit diagonal, respectively, 
and D is a diagonal matrix, is used for solving matrix equations AX = B. We construct 
a similar decomposition for the Bellman equation X = AX © B. 

For the case AX = B, the decomposition A = LDM induces the following decomposi- 
tion of the initial equation: 

(20) LZ = B, DY = Z, MX = Y. 
Hence, we have 

(21) A' 1 = M~ l D- l L- x , 

if A is invertible. In essence, it is sufficient to find the matrices L, D and M, since the 
linear system AX = B is easily solved by a combination of the forward substitution for 
Z, the trivial inversion of a diagonal matrix for Y, and the back substitution for X. 

Using the LDM-factorization of AX = B as a prototype, we can write 

(22) Z = LZ®B, Y = DY®Z, X = MX © Y. 
Then 

(23) A* = M*D*L*. 

A triple (L,D,M) consisting of a lower triangular, diagonal, and upper triangular 
matrices is called an LDM-factorization of a matrix A if relations (T221 and (T231 are 
satisfied. We note that in this case, the principal diagonals of L and M are zero. 

Our universal modification of the LDM-factorization used in matrix analysis for the 
equation AX = B is similar to the L?7-factorization of Bellman equation suggested by 
Carre in 

If A is a symmetric matrix over a semiring with a commutative multiplication, the 
amount of computations can be halved, since M and L are mapped into each other under 
transposition. 

We begin with the case of a triangular matrix A = L (or A = M). Then, finding 
X is reduced to the forward (or back) substitution. Note that in this case, equation 
X = AX © B has unique solution, which can be found by the obvious algorithms given 
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below. In these algorithms B is a vector (denoted by b), however they could be modified 
to the case when B is a matrix of any appropriate size. We are interested only in the 
case of strictly lower-triangular, resp. strictly upper-triangular matrices, when = for 
i < ji resp. a,ij = for i > j. 

Algorithm 5. Forward substitution. 

Input: Strictly lower-triangular n x n matrix I] 
n x 1 vector b. 

for k = 2 : n 

y{k) = l{k,l : k- 1) Q y{l : k- 1) 
end 

Output: vector y. 

Algorithm 6. Backward substitution. 

Input: Strictly upper-triangular n x n matrix m; 
fix 1 vector b. 

for k — n — 1 : —1 : 1 

y(k) = m(k, k + 1 : n) y{k + 1 : n) 
end 

Output: vector y. 

Both algorithms require n 2 /2 + 0(n) operations © and 0, and no algebraic closures. 

After performing a LDM-decomposition we also need to compute the closure of a diag- 
onal matrix: this is done entry wise. 

We now proceed with the algorithm of LDM decomposition itself, that is, computing 
matrices L, D and M satisfying (122]) and (1231) . First we give an algorithm, and then we 
proceed with its explanation. 

Algorithm 7. LDM-decomposition (version 1). 

Input: an n x n matrix A with entries a(i,j), 

also used to store the final result 

and intermediate results of the computation process. 

for j — 1 : n — 1 

a(j + 1 : n, j) = a(j + 1 : n, j) v(j) 

a(j + 1 : n,j + 1 : n) = a(j + 1 : n,j + 1 : n) © a(j + 1 : n,j) a(j,j + 1 : n) 
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a(j,j + 1 : n) = v(j) a(j,j + l:n) 



end 



The algorithm requires n 3 /3 + 0(n 2 ) operations © and 0, and n — 1 operations of algebraic 
closure. 

The strictly triangular matrix L is written in the lower triangle, the strictly upper 
triangular matrix M in the upper triangle, and the diagonal matrix D on the diagonal of 
the matrix computed by Algorithm [71 We now show that A* = M*D*L*. Our argument 
is close to that of [3]. 

We begin by representing, in analogy with the escalator method, 



as the multiplication on the r.h.s. leads to expressions fully asnalogous to (TIT]) , where 
(h^a^g^ plays the role of Uk+i- Here and in the sequel, Okxi denotes the k x I 

matrix consisting only of zeros, and denotes the identity matrix of size I. This can be 
also rewritten as 



(24) 




It can be verified that 




(25) 




(26) 



A* = M*Dl{A (2) y LI, 



where 
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/ o lxl 1 x(n— 1) 

1 — * (11 r\ 

\ a ll9 u (n-l)x(n-l) 



(27) =ftW fl ;/)©BW. 

Here we used in particular that L 2 = and M 2 = and hence L\ = I © Li and 
Aff = / ©Mi. 

The first step of Algorithm [7J (k = 1) computes 



(28) f ^eL.eM.ec,, 



which contains all relevant information. 

We can now continue with the submatrix of A^ 2 ^ factorizing it as in ( 13. 3 p and (126]) . 
and so on. Let us now formally describe the fcth step of this construction, corresponding 
to the kth step of Algorithm [71 On that general step we deal with 



(29) A® 



-l)x(fc-l) -l)x(n— fc+1) 

0(n-fc+l)x(Jfc-l) 



where 



4? ^ (fc) 

^(*) #(*) 



(30) 

Like on the first step we represent 
(31) = M^D^A^ 1 ))*! 



where 
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(32) 



D, 



R (k+i) 



(O 



(fc-l)x(fc-l) 

Olx(fc-l) 



0(fe-l)xl 

o lxl 



o 



(fc-l)x(n-fc) 



\0( n _fc) x (fc-l) 0(n-fe)xl 0(n-k)x{n-k)J 

l)x(fc-l) 0(fc_l) x l 0(fc- 

(fc-1) a fcfc (n-fe) 

\0( n _fc) x (fc-l) C( n _fe)xi 0(„_fc) x (n-fe)/ 

l)x(fc-l) 0(fc_l) x l 

01 x (fc-1) Oixl Olx(n-fc) 

\0(n-fc)x(fc-l) {a>k k )*9^ 0(n-k)x(n-k)J 

Okxk Ofcx(n— fc) 
0( n -k)xk R^ k+1 ^ 



Note that we have the following recursion for the entries of A^ h ': 



(33) 



ffc+i) 



0. 



if i < k or j < k, 



o^eog^agj)^, otherwise. 



This recursion is immediately seen in Algorithm [71 Moreover it can be shown by induction 
that the matrix computed on the kth step of that algorithm equals 



(34) 



A (*+i) e L . e M . e A . 



i=l 



i=l 



8=1 



In other words, this matrix is composed from h^a^, h^ k \a kk )* (in the upper triangle), 

^on the diagonal), and R^ k+1 ^ 



(k) 



' a kk 



a n g w , (a kk )*g^ > (in the lower triangle), an, 
(in the south-eastern corner). 

After assembling and unfolding all expressions (13TT) for A {k \ where k = 1, . . . , n, we 
obtain 



(35) 



A* = M*D{ ■ ■ ■ M*D* n L* n ■■■L\. 



(actually, M n = L n = and hence M* = L* = I). Noticing that D* and M* commute 
for i < j we can rewrite 



(36) 



A* = M* - • • M*Dl ■ ■ ■ D* n L* n ■•■L* v 
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Consider the identities 

(£>!©...© D n y = £>*••■ D* n , 
(37) (L 1 ®...®L n )* = L* n ---L{, 

(M 1 @...@M n y = M*---M*. 

The first of these identities is evident. For the other two, observe that M| = Li = for 
all k, hence M£ / M k and /. A . / /./,. Further, I^Zy = for % > j and = 

for i < j. Using these identities it can be shown that 

re-l 

(Li © . . . © L n )* = 0(Li © ... © L n ) i = 

i=0 

= (/©L n )---(/ffiL 1 ) = L;---Lt, 

(38) 

v 7 n-1 

(Mi © ... © M n )* = 0(Mi © ... © M n )* = 

i=0 

= (J © Mi) •••(/© M n ) = M* • ■ ■ M*, 

which yields the last two identities of ([3"T]) . Notice that in (138]) we have used the nilpotency 
of L x © . . . © L n and Mi © ... © M n , which allows to apply (jSJ). 

It can be seen that the matrices M := Mi © ... © M n , L := L x © . . . © L n and D : = 
Dx © . . . © _D n are contained in the upper triangle, in the lower triangle and, respectively, 
on the diagonal of the matrix computed by Algorithm [71 These matrices satisfy the LDM 
decomposition A* = M*D*L*. This concludes the explanation of Algorithm [7J 

In terms of matrix computations, Algorithm [7J is a version of LDM decomposition with 
outer product. This algorithm can be reorganized to make it almost identical with |14j . 
Algorithm 4.1.1: 

Algorithm 8. LDM- decomposition (version 2). 

Input: an n x n matrix A with entries a(i,j), 

also used to store the final result 

and intermediate results of the computation process. 

for j — 1 : n 

v(l : j) = a(l : j,j) 

for k = 1 : j - 1 v (k + 1 : j) = v(k + 1 : j) © a(k + 1 : j, k) v(k) 
end 

for i — 1 : j — 1 

a(i,j) = (a(i,i))* Q v(i) 
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end 



a(j,j) = v{j) 
for k = 1 : j — 1 

a(j + 1 : to, j) = a(j + 1 : to, j) © a(j + 1 : to, k) 
end 

d = Mi))* 

a(j + 1 : to, j) = a(j + 1 : to, j) d 



This algorithm performs exactly the same operations as Algorithm [71 computing consec- 
utively one column of the result after another. Namely, in the first half of the main loop 
it computes the entries ay for i = 1, . . . ,j, first under the guise of the entries of v and 
finally in the assignment u a(i,j) = (a(i,i))* © v (i)". In the second half of the main loop 
it computes a$ . The complexity of this algorithm is the same as that of Algorithm [7J 

3.4. LDM decomposition with symmetry and band structure. When matrix A is 
symmetric, i.e., = aji for all it is natural to expect that LDM decomposition must 
be symmetric too, that is, M = L T . Indeed, going through the reasoning of the previous 
section, it can be shown by induction that all intermediate matrices A^ k > are symmetric, 
hence M k = L T k for all k and M = L T . We now present two versions of symmetric 
LDM decomposition, corresponding to the two versions of LDM decomposition given in 
the previous section. Notice that the amount of computations in these algorithms is 
nearly halved with respect to their full versions. In both cases they require to 3 /6 + 0(n 2 ) 
operations © and ©(each) and to — 1 operations of taking algebraic closure. 

Algorithm 9. Symmetric LDM- decomposition (version 1). 

Input: an to x to symmetric matrix A with entries a(i,j), 

also used to store the final result 

and intermediate results of the computation process. 

for j — 1 : to — 1 



end 



- j + 1 : to 
j + l:k 



for k 
for I -- 

a(k, I) 



a{k, 1) © a{k, j) v(j) a(l, j) 



end 



end 



a(j + 1 : n,j) 
end 
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= a(j + 1 : n,j) v(j) 
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The strictly triangular matrix L is contained in the lower triangle of the result, and the 
matrix D is on the diagonal. 

The next version generalizes [2] Algorithm 4.1.2. Like in the prototype, the idea is to 
use the symmetry of A precomputing the first j — 1 entries of v inverting the assignment 
"a(i, j) = a(i, i)* v(i)" for i — 1, . . . ,j — 1. This is possible since a(j, i) = a(i,j) belong 
to the first j — 1 columns of the result that have been computed on the previous stages. 

Algorithm 10. Symmetric LDM- decomposition 
(version 2). 

A is an n x n symmetric matrix with entries a(i,j), 

also used to store the final result 

and intermediate results of the computation process. 

for j — 1 : n 
for i — 1 : j — 1 

v{i) = (a{i ) i)*y l a{i ) i) 
end 

v (j) = v (j) © a {j, 1 : j - 1) v(l : j - 1) 
a(j,j) = v(j) 
for k = 1 : j — 1 

a(j + 1 : n, j) = a(j + 1 : n, j) © a(j + 1 : n, k) ?;(&) 
end 

d = (v(j)Y 

a(j + 1 : n, j) = a(j + 1 : n, j) Q d end 

Note that this version requires invertibility of the closures a(i,i)* computed by the algo- 
rithm. 

Remark 3. In the case of idempotent semiring we have (D*) 2 = D* , hence 
A* = (M* D*)(D*L*). When A is symmetric we can write A* = (G*) T G* where G = 
D*L. Evidently, this idempotent Cholesky factorization can be computed by minor 
modifications of Algorithms [PI and ITU See also [H] : Algorithm 4-2.2. 

A = (aij) is called a band matrix with upper bandwidth q and lower bandwidth p 
if aij = for all j > i + q and all i > j + p. A band matrix with p = q = 1 is called 
tridiagonal. To generalize a specific LDM decomposition with band matrices, we need 



30 G. L. LITVINOV, A. YA. RODIONOV, S. N. SERGEEV, AND A. N. SOBOLEVSKI 

to show that the band parameters of the matrices , . . . , A^ n > computed in the process 
of LDM decomposition are not greater than the parameters of A^ = A. Assume by 
induction that A = A^\ . . . , A^ have the required band parameters, and consider an 
entry a\- for % > j + p. If i < k or j < k then a\- = 0, so we can assume i > k and 

(k) 

j > k. In this case i > k + p, hence a ik = and 

(fc+l) (k) _ {k) f {k)s* (k) n 

Thus we have shown that the lower bandwidth of A^ is not greater than p. It can be 
shown analogously that its upper bandwidth does not exceed q. We use this to construct 
the following band version of LDM decomposition, see [H] Algorithm 4.3.1 for a prototype. 

Algorithm 11. LDM decomposition of a band matrix. 

A is an n x n band matrix with entries a(i,j), 

lower bandwidth p and upper bandwidth q 

also used to store the final result 

and intermediate results of the computation process. 

for j — 1 : n — 1 

v(j) = (a(j,j))* 

for % = j + 1 : min(j + p, n) 

a(ij) = a(i,j)®v(j) 

end 

for k = j + 1 : min(j + q, n) 
for i — j + 1 : min(j + p, n) 
a (k, j) = a(k, j) © a(k, i) a(i, j) 
end 
end 

for k = j + 1 : min(j + q, n) 
a(j, k) = v(j)Qa(j, k) 
end 
end 

When p and q are fixed and n » p, q is variable, it can be seen that the algorithm 
performs approximately npq operations © and © each. 

Remark 4. There are important special kinds of band matrices, e.g., Hessenberg and 
tridiagonal matrices. Hessenberg matrices are defined as band matrices with p = 1 and 
q = n, while in the case of tridiagonal matrices p — q — 1. It is straightforward to write 
further adaptations of Algorithm [TT] to these cases. 
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3.5. Iteration schemes. We are not aware of any truly universal scheme, since the 
decision when such schemes work and when they should be stopped depends both on the 
semiring and on the representation of data. 

Our first scheme is derived from the following iteration process: 

(39) X (fc+1) = AX {k) © B 

trying to solve the Bellman equation X = AX © B. Iterating expressions f !39p for all k 
up to m we obtain 

TO— 1 

(40) X (m) = A m X {0) © A*B 

i=Q 

Thus the result crucially depends on the behaviour of A m X^\ The algorithm can be 
written as follows (for the case when B is a column vector). 

Algorithm 12. Jacobi iterations 

Input: n x n matrix A with entries a(i,j); 
n x 1 column vectors b and x 

situation= 'proceed' 

while situation=='proceed' 

x = A x © b 

situation=newsituation(. . . ) 

if situation=='no convergence' 

disp('Jacobi iterations did not converge') 

exit 

end 

if situation=='convergence' 

disp(' Jacobi iterations converged') 

exit 

end 

end 

Output: situation, x. 

Next we briefly discuss the behaviour of Jacobi iteration scheme over the usual arithmetic 
with nonnegative real numbers, and over semiring R max . For simplicity, in both cases we 
restrict to the case of irreducible matrix A, i.e., when the associated digraph is strongly 
connected. 
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Over the usual arithmetic, it is well known that (in the irreducible nonnegative case) 
the Jacobi iterations converge if and only if the greatest eigenvalue of A, denoted by r(A), 
is strictly less than 1. This follows from the behaviour of A m x^°\ In general we cannot 
obtain exact solution of x = Ax + b by means of Jacobi iterations. 

In the case of R m ax, the situation is determined by the behaviour of A m x^ which 
differs from the case of the usual nonnegative algebra. However, this behaviour can be 
also analysed in terms of r(A), the greatest eigenvalue in terms of max-plus algebra (that 
is, with respect to the max-plus eigenproblem A x = A x). Namely, A m x^ — > and 
hence the iterations converge if r(A) < 1. Moreover A* = (I © A © ... © A n ~ l ) and hence 
the iterations yield exact solution to Bellman equation after a finite number of steps. 
To the contrary, A m x^ — > +oo and hence the iterations diverge if r(A) > 1. See e.g. [7] 
for more details. On the boundary r(A) = 1, the powers A m reach a periodic regime after 
a finite number of steps. Hence A*b © A m x^ also becomes periodic, in general. If the 
period of A m x^ is one, i.e., if this sequence stabilizes, then the method converges to a 
general solution of x = Ax © b described as a superposition of A*b and an eigenvector of 
A P, [22]. The vector A*b may dominate, in which case the method converges to A*b as 
"expected" . However, the period of A*b © A m x^ may be more than one, in which case 
the Jacobi iterations do not yield any solution of x = Ax © b. See [5] for more information 
on the behaviour of max-plus matrix powers and the max-plus spectral theory. 

In a more elaborate sheme of Gauss-Seidel iterations we can also use the previously 
found coordinates of X^ k \ In this case matrix A is written as L © U where L is the 
strictly lower triangular part of A, and U is the upper triangular part with the diagonal. 
The iterations are written as 

(41) X (fe) = LX (fc) © UX {k - l) ®B = L*UX^ 1] © L*B 

Note that the transformation on the r.h.s. is unambiguous since L is strictly lower trian- 
gular and L* is uniquely defined as / © L © ... © L n_1 (where n is the dimension of A). 
In other words, we just apply the forward substitution. Iterating expressions (jUJ) for all 
k up to m we obtain 

171 — 1 

(42) X (m) = {L*U) m X {0) © (£)(L*UyL*B 
The r.h.s. reminds of the formula [L © U)* = 

(L*U)*L*, see ([Tj), so it is natural to expect that these iterations converge to A*B with a 
good choice of X^°\ The result crucially depends on the behaviour of (L*U) m X^ . The 
algorithm can be written as follows (we assume again that B is a column vector). 
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Algorithm 13. Gauss-Seidel iterations 
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Input: n x n matrix A with entries a(i,j); 
n x 1 column vectors b and x 

situation= 'proceed' 

while situation=='proceed' 

for i = 1 : n 

y(i) = a(i, i : n) x(i : n) © 
end 

for % = 2 : n 

x(i) = a(i, 1 : i — 1) x(l : i — 1) 
' end 

situation= newsituation( . . . ) 

if situation=='no convergence' 

disp('Gauss-Seidel iterations did not converge') 

exit 

end 

if situation=='convergence' 

disp('Gauss-Seidel iterations converged') 

exit 

end 

end 

Output: situation, x. 

It is plausible to expect that the behaviour of Gauss-Seidel scheme in the case of max-plus 
algebra and nonnegative linear algebra is analogous to the case of Jacobi iterations. 

3.6. Software implementation of universal algorithms. Software implementations 
for universal semiring algorithms cannot be as efficient as hardware ones (with respect to 
the computation speed) but they are much more flexible. Program modules can deal with 
abstract (and variable) operations and data types. Concrete values for these operations 
and data types can be defined by the corresponding input data. In this case concrete 
operations and data types are generated by means of additional program modules. For 
programs written in this manner it is convenient to use special techniques of the so- 
called object oriented (and functional) design, see, e.g., [35j H3j H9] . Fortunately, powerful 
tools supporting the object-oriented software design have recently appeared including 
compilers for real and convenient programming languages (e.g. C + + and Java) and 
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modern computer algebra systems. Recently, this type of programming technique has 
been dubbed generic programming (see, e.g., [43J). 

C + + implementation Using templates and objective oriented programming, Churkin 
and Sergeev [50] created a Visual C + + application demonstrating how the universal 
algorithms calculate matrix closures A* and solve Bellman equations x = Ax © b in 
various semirings. The program can also compute the usual system Ax = b in the usual 
arithmetic by transforming it to the "Bellman" form. Before pressing "Solve", the user 
has to choose a semiring, a problem and an algorithm to use. Then the initial data are 
written into the matrix (for the sake of visualization the dimension of a matrix is no more 
than 10 x 10). The result may appear as a matrix or as a vector depending on the problem 
to solve. The object-oriented approach allows to implement various semirings as objects 
with various definitions of basic operations, while keeping the algorithm code unique and 
concise. 

Examples of the semirings. The choice of semiring determines the object used by the 
algorithm, i.e., the concrete realization of that algorithm. The following semirings have 
been realized: 

1) © = + and © = x: the usual arithmetic over reals; 

2) © = max and © = +: max-plus arithmetic over RU {— oo}; 

3) © = min and © = +: min-plus arithmetic over R U {+oo}; 

4) © = max and © = x : max-times arithmetic over nonnegative numbers; 

5) © = max and © = min: max-min arithmetic over a real interval [a, b] (the ends a 
and b can be chosen by the user); 

6) © =OR and © =AND: Boolean logic over the two-element set {0, 1}. 

Algorithms. The user can select the following basic methods: 

1) Gaussian elimination scheme, including the universal realizations of escalator 
method (Algorithm [1]), Floyd- Warshall (Algorithm [21 Yershov's algorithm (based 
on a prototype from [TTJ Ch. 2), and the universal algorithm of Rote [44] ; 

2) Methods for Toeplitz systems including the universal realizations of Durbin's 
and Levinson's schemes (Algorithms [3] and HJ) ; 

3. LDM decomposition (Algorithm [7j) and its adaptations to the symmetric case 
(Algorithm [9j, band matrices (Algorithm [TT]) . Hessenberg and tridiagonal matri- 
ces. 

4) Iteration schemes of Jacobi and Gauss-Seidel. As mentioned above, these 
schemes are not truly universal since the stopping criterion is different for the 
usual arithmetics and idempotent semirings. 
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Types of matrices. The user may choose to work with general matrices, or with a matrix 
of special structure, e.g., symmetric, symmetric Toeplitz, band, Hessenberg or tridiagonal. 

Visualization. In the case of idempotent semiring, the matrix can be visualized as a 
weighted digraph. After performing the calculations, the user may wish to find an optimal 
path between a given pair of nodes, or to display an optimal paths tree. These problems 
can be solved using parental links like in the case of the classical Floyd- Warshall method 
computing all optimal paths, see, e.g., |45J. In our case, the mechanism of parental links 
can be implemented directly in the class describing an idempotent arithmetic. 

Other arithmetics and interval extensions. It is also possible to realize various types 
of arithmetics as data types and combine this with the semiring selection. Moreover, all 
implemented semirings can be extended to their interval versions. Such possibilities were 
not realized in the program of Churkin and Sergeev [50], being postponed to the next 
version. The list of such arithmetics includes integers, and fractional arithmetics with the 
use of chain fractions and controlled precision. 

MATLAB realization. The whole work (except for visualization tools) has been du- 
plicated in MATLAB [50], which also allows for a kind of object-oriented programming. 
Obviously, the universal algorithms written in MATLAB are very close to those described 
in the present paper. 

Future prospects. High-level tools, such as STL [13J 09j, possess both obvious advan- 
tages and some disadvantages and must be used with caution. It seems that it is natural 
to obtain an implementation of the correspondence principle approach to scientific cal- 
culations in the form of a powerful software system based on a collection of universal 
algorithms. This approach should ensure a working time reduction for programmers and 
users because of the software unification. The arbitrary necessary accuracy and safety of 
numeric calculations can be ensured as well. 

The system has to contain several levels (including programmer and user levels) and 
many modules. 

Roughly speaking, it must be divided into three parts. The first part contains modules 
that implement domain modules (finite representations of basic mathematical objects). 
The second part implements universal (invariant) calculation methods. The third part 
contains modules implementing model dependent algorithms. These modules may be used 
in user programs written in C ++ , Java, Maple, Matlab etc. 

The system has to contain the following modules: 

— Domain modules: 

— infinite precision integers; 
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— rational numbers; 

— finite precision rational numbers (see [16]); 

— finite precision complex rational numbers; 

— fixed- and floating-slash rational numbers; 

— complex rational numbers; 

— arbitrary precision floating-point real numbers; 

— arbitrary precision complex numbers; 

— p-adic numbers; 

— interval numbers; 

— ring of polynomials over different rings; 

— idempotent semirings; 

— interval idempotent semirings; 

— and others. 
— Algorithms: 

— linear algebra; 

— numerical integration; 

— roots of polynomials; 

— spline interpolations and approximations; 

— rational and polynomial interpolations and approximations; 

— special functions calculation; 

— differential equations; 

— optimization and optimal control; 

— idempotent functional analysis; 

— and others. 

This software system may be especially useful for designers of algorithms, software 
engineers, students and mathematicians. 
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